 %%
str = 'noise0.01.bmp';
global HD
global HR
HD = frameFilterD('linear');
HR = frameFilterR('linear');
I = readImage( str );
% I = imnoise(I,'gaussian',0,0.01);
%%
nPhases = 4;
nu = 1;
beta = 1.5;
[c label] = Kmeans(I, nPhases, 50);
c = sort(c);
M = size(I,1);
N = size(I,2);
f = zeros(M,N,nPhases);
for i = 1 : nPhases
    f(:,:,i) = nu*(abs(I - c(i))).^beta;
end
nLevels = 2;
u = zeros(M,N,nPhases);
for i = 1 : nPhases;
    y{i} = frameDec2(I,nLevels,HD);
%     u(:,:,i) = (label == i);
end
ubar = u;
%%
nIterations = 10;
tic;
delta = 0.2;
tau = 0.2;
for iter = 1 : nIterations
    for i = 1 : nPhases
        tmp = frameDec2(ubar(:,:,i), nLevels, HD);
        y{i} = cellPlus(y{i}, cellMultiply(tmp, delta));
        y{i} = thresholdCell_2(y{i},1,HD);
        u(:,:,i) = u(:,:,i) - tau*frameRec2(y{i}, nLevels, 0, HR) -tau*f(:,:,i);
    end
    ubar = u;
    u = shiftdim(projSplx(shiftdim(u,2)),1);
end
[Y U] = max(u,[],3);
output = c(U);

figure;
% imagesc(I,[0,1]); axis off; axis equal;
% colormap(gray);hold on;
% contour(U,[4,3],'r');
toc
 imshow(output,[])




















